pacman::p_load(dplyr, ggplot2, vcd, plotly)

【A】實務案例一：擲筊的成功機率

🧙 問題討論：
在cup.csv檔案裡面是一千個信徒的擲筊次數，假定每一個人都是擲到3次成功才停止，請問：
  ■ 這些擲筊次數是甚麼分佈？ 它的參數是？
  ■ 這個筊的成功機率大約是？
  ■ 請畫出用這個筊擲10次成功之前，失敗次數的機率分布？
  ■ 用這個筊擲15次還不能有5次成功的機率是？

失敗次數是Poisson嗎？ (因為擲到3次成功才停止，所以失敗次數等於擲筊次數減三)

X = read.csv("data/cup.csv")$x
fit =  goodfit(table(X - 3), "poisson") #NO
summary(fit)

     Goodness-of-fit test for poisson distribution

                    X^2 df    P(> X^2)
Likelihood Ratio 1605.4 33 5.0416e-317

已知成功次數為3，用Negative Binomial試試看

fit =  goodfit(table(X - 3), "nbinom", par=list(size=3)) #YES
summary(fit)

     Goodness-of-fit test for nbinomial distribution

                    X^2 df P(> X^2)
Likelihood Ratio 37.564 33  0.26795
fit$par  #par=list(size=3)
$size
[1] 3

$prob
[1] 0.248
         #參數，已知為3所以直接設定

🗿 : 這個筊的成功機率大約是？

fit$par$prob  #0.248
[1] 0.248

🗿 : NBinomial[n=3,p=0.248]是 … 的分佈？

#筊的成功機率大約是0.248的情況下，累積成功3次前所需丟擲次數(失敗次數)

#幾何分布在R定為首次成功之前所經歷的失敗次數

dnbinom(0:45, 3 , fit$par$prob) %>% barplot(names=0:45)

🗿 : 畫出用這個筊擲10次成功之前，失敗次數的機率分布

par=c(margin=c(3,3,3,1),cex=0.7)
dnbinom(0:80, 10 , fit$par$prob) %>% barplot(names=0:80)

🗿 : 用這個筊擲15次還不能有5次成功的機率是？

1 - pnbinom(10,5,fit$par$prob)
[1] 0.69309
#丟15次，5次成功前失敗次數大於10的機率

pbinom(4,15,fit$par$prob)
[1] 0.69309
#成功機率0.248時，成功次數小於等於4

【B】實務案例二：老忠實賭局

🗿 問題討論：
假如賭場老闆從零開始把每10秒鐘設為一個區間，每個區間的賭金是五塊錢 …
 ■ 你要怎麼押注，才能獲得最高的期望報酬呢？
 ■ 你的賭金和期望報酬各是多少？

§ 使用R的內建的功能建立PDF模型

D = faithful$eruptions        # copy data to a short name
Adjust = 0.5                       # bandwidth adjustment
DEN = density(D, adjust = Adjust)      # density function
PDF = approxfun(DEN$x, DEN$y, yleft=0, yright=0)    # PDF

par(mar=c(2,4,2,1), cex=0.7)
curve(PDF(x), 1, 5.6, col='blue',lwd=2,main="Prob. Density Function")
rug(D)
abline(h=seq(0,1,0.1),col='lightgrey',lty=3)

§ 利用模型策略規劃

x1 = seq(1, 6, 1/6); x2 = x1 + (1/6)     # 臨界值
px = sapply(1:length(x1), function(i) {  # 機率
  integrate(PDF, x1[i], x2[i])$value     #       
  })  
net = 100 * px - 5                       # 淨期望報酬
df = data.frame(x1, x2, px, net) %>% 
  arrange(desc(net)) %>%
  mutate(
    payoff =  cumsum(net),               # 累計(淨)期望報酬
    investment = row_number() * 5,       # 累計投資、資本需求
    roi = payoff / investment            # 累計投報率
    ) 
df %>% filter(net > 0) %>% round(3)
     x1    x2    px   net payoff investment   roi
1 4.333 4.500 0.096 4.577  4.577          5 0.915
2 4.500 4.667 0.089 3.938  8.515         10 0.852
3 4.167 4.333 0.089 3.932 12.447         15 0.830
4 1.833 2.000 0.082 3.226 15.673         20 0.784
5 4.000 4.167 0.077 2.654 18.327         25 0.733
6 2.000 2.167 0.070 1.968 20.294         30 0.676
7 4.667 4.833 0.068 1.843 22.137         35 0.632
8 1.667 1.833 0.064 1.367 23.504         40 0.588
9 3.833 4.000 0.058 0.834 24.338         45 0.541
nrow(df %>% filter(net > 0))         #要押幾注，Ans:9注，前提要有45塊
[1] 9
5*nrow(df %>% filter(net > 0))       #投資金額
[1] 45
sum(filter(df,net > 0)$payoff)       #最高期望報償
[1] 149.81

🗿 問題討論：
將獲利的期望報酬和賭金的比值稱為「期望投資報酬率」 …
 ■ 「最大期望投資報酬率」的策略是？ A:4.333 ~ 4.500
 ■ 這個策略的期望報酬和報酬率各是多少？ A:(淨)期望報酬:4.577 & ROI:0.915
 ■ 它和最高期望報酬的策略是相同的嗎？ A:不是
 ■ 哪一個策略目標比較好呢？ A:各有千秋
 ■ 我們應該如何選擇(設定)策略目標呢？ A:至少選擇淨期望報酬>0的目標 

df %>% mutate(strategy = case_when(
  payoff < 0 ~ "虧損策略",
  net < 0 ~ "無效率策略",
  TRUE ~ "有效率策略"
  )) %>% 
  ggplot(aes(x=payoff, y=roi,col=strategy)) + 
  geom_point(size=3) +
  xlab("期望獲利") + ylab("投資報酬率") + ggtitle("策略空間") +
  theme_bw() -> g
ggplotly(g)

💡 學習重點：
  ■ 如何求出機率分布函數？
    A. 將實證資料fit進已知的理論分布、並求出參數
    B. 直接就實證資料做成(平滑)機率密度函數
  ■ 即使是很簡單的一個案例，我們也需要綜合使用機率、程式和商業知識，才能做出合適的決策！
  ■ 機率分布的商業應用套路：
    1. 寫出各事件(變數的值或區間)的報償
    2. 求出變數的機率分布函數
    3. 求出各事件的機率
    4. 求出不同策略的期望報償
    5. 依公司的策略目標做選擇